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Abstract 

5jry This paper concerns a fast, one-step iterative technique of imaging extended perfectly con- 
ducting cracks with Dirichlet boundary condition. In order to reconstruct the shape of cracks 
from scattered field data measured at the boundary, we introduce a topological derivative- 
based electromagnetic imaging function operated at several nonzero frequencies. The proper- 
ties of the imaging function are carefully analyzed for the configurations of both symmetric 
and non-symmetric incident field directions. This analysis explains why the application 
of incident fields with symmetric direction operated at multiple frequencies guarantees a 
successful reconstruction. Various numerical simulations with noise-corrupted data are con- 
ducted to assess the performance, effectiveness, robustness, and limitations of the proposed 
technique. 
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1. Introduction 

One goal of the inverse scattering problem is to find the shape of bulk or flaw in a 
medium using the scattered field data measured at the boundary. This is an old but it 
has several of applications such as detection of cracks in material engineering, ultrasound 
imaging in medical sciences, and scanning anti-personnel mines hidden in the ground in 
military services. Related researches may be found in 0, Q, [l6, 17, 18, 37] and the references 



therein. In many studies, various shape reconstruction algorithms have been developed, 
most of which are focused on the minimization of the least-square functional by using a 
Newton-type iteration strategy, e.g., level-set method 0, [l8, 35]. The merits of the iteration 



strategy is that it does not require a very large number of boundary measurements for 
complete shape reconstruction; instead, in order to guarantee successful reconstruction, it 
demands very high computational expenditures, optimal regularization terms related to the 
problems at hand, calculation of a complex Frechet derivative at each iteration, and a priori 
information of the unknown target to be reconstructed. Nevertheless, even if these conditions 
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are satisfied, the iteration strategy must begin with an initial guess that is close to the true 
one in order to avoid undesirable situations such as non-convergence or falling into the local 
minima. Therefore, development of a mathematical theory and an algorithm for generating 
a good initial guess is an important research topic. 

Conversely, certain non-iterative shape reconstruction algorithms have been proposed; for 
example. Multiple Signal Classification (MUSIC) algorithm 0, E3, 32, 34 ] , linear sampling 
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example. Multiple signal uiassmcation (lyiusiu) aigontnm p, LUJ, [dz|, lojj , linear samp 
method 15. 21 . Kirchhoff migration 0, [28, 29, 33 , and Fourier inversion based one 0,0, 
In contrast with the iterative strategy, these algorithms require a large amount of incident 
field data and a significant number of boundary measurements. However, if these conditions 
are satisfied, non-iterative shape reconstruction algorithms can prove to be effective, fast, 
and easy to extend to multiple targets. 

Topological derivative concept is a non-iterative strategy. This concept was originally 
developed for the shape optimization problem, but its application to rapid shape reconstruc- 



tion has only recently been proven. Related works can be found in [5|, [10|, [T2|, |13j, |14j, |24 



25l . 1271 . l30l l3ll |38| and references therein. One of the advantages of topological derivative 



concept is that it does not require a large amount of many incident field data; however, a 
reduction in the amount of this data has been reported to result in poor resolution of the 
reconstructed shape. 

The purpose of this literature is to establish an effective and fast reconstruction algorithm 
for detecting both the location and shape of perfectly conducting crack(s) from boundary- 
measured scattered field data by using only a small amount of incident field data. Based 
on the structure of topological derivative imaging functional operated at a fixed-frequency 
24| and multi-frequency imaging schemes (fl |26| , we introduce a multi-frequency topological 
derivative concept. Note that recent work |4|, Formula (1.33)] has briefly introduced a multi- 
frequency topological derivative but through a more detailed analysis, we investigate the 
structure of the proposed multi-frequency topological derivative and its benefits over the 
traditional topological derivative. 

The remainder of this paper is organized as follows. Section |2] briefly introduces the 
traditional topological derivative by creating a linear, narrow crack. Section |3] describes 
the design of a multi-frequency topological derivative imaging function and analyzes it to 
investigate its properties. Section H] presents various numerical experiments to assess the 
performance of the proposed imaging function under various situations. Section contains 
a short conclusion and some remarks on future work. 



2. Topological derivative: an inspection 

Assume that a perfectly conducting crack C is completely hidden in a homogeneous 
domain DcK 2 with a smooth boundary B. We denote x c as a two-dimensional vector that 
lies on crack C. 

Throughout this paper, we consider the so-called Transverse Magnetic (TM) polarization 
case, letting w^ tal (x; k) be the (single-component) electric field that satisfies the following 
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boundary value problem: 

' Au^ tal (£; k) + k 2 u^\x; k) = in V\C 



u 



total . 



x;k) = on C 



du. 



total i 



X] k) dexp(ik6 m ■ x) 



on B, 



1, 2, • • • , M, denotes 



dv(x) dv(x) 

where v(x) represents the unit outward normal to x G B and 9 m , m 
the two-dimensional vector on the unit circle S 1 . Throughout this paper, we let k denote a 
strictly positive wavenumber and presume k 2 to not be an eigenvalue of ([1]). Similarly, let 
x] k) = exp(ik8 m ■ x) be the solution of equation ([!]) without C. Then, the problem 



u 



back. 



we consider here is the computation of the topological derivative of the residual function 
depending on the solution u^ t3l (x; k): 

M M 



(V-k) -.= \^\\uT\^k)-v^ c \x-k) 



m=l 



m=l 



U 



total / 



(£; k) -< ack (f; fc)| dB(x), 

(2) 

where |it^ tal (f; k)\ 2 = u^ ta \x; k)u^ tal (x; k). 

In order to compute the topological derivative, let us create a small linear crack C of 
length 2£(<^C k) centered at x s = (xl,xf) G V\B such that 

£ = {(t 1 x 2 s ):x 1 s -e<z<x 1 s + e}, 

and denote T> V£ as that domain. Then, because of the change in the topology of T>, we can 
consider the topological derivative Ktd(^s! k) based on the residual function R(D; k) with 
respect to point follows: 

i(V V £; k) - R(X>; k) 



^TD 



x s ] k) 



lim 



(3) 



m k) 

where <j>{l] k) — > as £ — > 0+. From relationship (J3]), we have the following asymptotic 
expansion: 

R(V V £; k) = R(V; k) + (/>(£; A;)R TD (£ S ; k) + o{(j){£] k)). 

Then, the topological derivative Rtd(^; k) with M different incident waves at a given 
wavenumber k is as follows (see for its derivation). 

Lemma 2.1 (Topological derivative). Let 3?(/) denote the real part of f. Then, the 
topological derivative corresponding to (TJ|) is given by 

' M 

x s ; k) — 3? ^ ] v^ nt (x s ; k)v}^ ck (x s ; k) , (4) 



m=l 



where vffl nt (-; k) satisfies the following adjoint problem: 

Av^(x; k) + k 2 v^ t (x;k) 
<9< djnt (f; k) 



in V 



du(x) 



u 



total / 



(x; k) - u^(x] k) on B 



(5) 



3. Multi-frequency topological derivative: introduction and analysis 



In this section, we introduce a multi-frequency topological derivative imaging function 
and explore its properties. Based on recent research 0,[lO, 30, 31], a sufficiently large number 
of incident directions are required to obtain a good result. This means that if M is small, an 
image with poor resolution will be reconstructed. Moreover, when boundary measurement 
data are affected by a large amount of random noise, (pQ) will fail to reconstruct a good image. 
Inspired by multi-frequency imaging algorithms (a. liil. 29 . |33j, we consider the following 
multi-frequency topological derivative for several wavenumbers {kf : f — 1, 2, • • • , F}: 



F ^ max |M TD (x s ; k f )\ 

J x a £V 



(6) 



where Mtd(^; kf) is given by for k — kf, f — 1, 2, • • • , F. 

In this section, we investigate the structure of ()6]). For this purpose, we recall a useful 
result (see jiij) as follows. 

Lemma 3.1. Let A ~ B imply that there exists a constant C such that A = BC . Then, 
for x s G V and x c G C, ^ satisfies 



'td(^ s ; k) ~ K 



M 



^ exp(ik6 m • (x c - x s ))$(x c , x s \ k) 



m=l 



(7) 



Here Q(x, y\ k) denotes the two-dimensional time-harmonic fundamental solution (or Green's 
function) of Helmholtz equation 



J {k\x-y\)+iY {k\x-y\)), 



<$>(x,y; fc) = --Hl{k\x-y\) 



where Hq, Jq, and Yq denote the Hankel function of the first kind, Bessel function, and 
Neumann function of order zero, respectively. 

Consequently, ([HD is analyzed for the following three cases of interest. 

Theorem 3.2 (Case 1: Symmetric incident directions). Under the configuration of a 
small number of symmetric incident directions, the structure of is 

M 



m=l 



X c - Xc 



\x c x s 



for a sufficiently large F and kp — > oo. 
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Proof. Assume that the symmetric incident directions are of the form: 

2(m-l)n . 2(m-l)n 



Or, 



COS 



-, sin ■ 



m = 1,2,- •• , M, 



M M 

where M is an even number, say, M = 2N. Based on this form, Q\ satisfies Q\ = —0i+n for 



1 = 1,2,--- ,N. 



In contrast with Lemma l3TT| we assume that M is not sufficiently large. Let ^s(f) denote 
the imaginary part of / then since 



3 



exp(ik9i ■ (x c - x s )) 
exp(ik9i ■ (x c - x s )) 



3? 



-9 



exp(ik9i +N ■ (x c - x s )) 
exp(ik0i +N ■ (x c -x s )) 



becomes 
M T d(x s ; k f ) ~ U 
= SR 



M 



— ^2 exp(%(9 m ■ (x c - x g ))$(x c , x s ; kf) 



J m=l 
N 



J- cos(k f 9 m ■ (x c - x s )){J (kf\x c -x s \) + iY (k f \x c - x s \)} 
2k f 

1 N 

— 22 cos(kfd m ■ (x c - x s ))J (kf\x c - x s \). 



2k f 



m=l 



Therefore, 



^td{x s ; kf 



max\R TD (x s ; k f )\ 



1 N 

— ^2 c os(k f 9 m ■ (x c - x s ))J (kf\x c - x s \) 



Let us recall that for > -1 (see Q, formula 11.4.37 (p. 487)]) 

1 



cos at J^(ht)dt = < 



Vb' 



|2 /~r2 



cos I fi sin — 



if < b < a 



b^ ( {in 

-. . — cos ( — ) if < a < b. 

^a 2 -b 2 (a + Va 2 - b 2 Y \ 2 



Let a = 9 m ■ (x c — x s ), b = \x c — x s \, and /i = 0. Then since 



h 2 2 — 1 "* ~* 1 2 




)•] 


>o, 
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by letting kp — > oo, (EJ) becomes 



mvitd(£. 



^TD 



%\ kf) 



N 



F max |M TD | N ^ J 

1 N 



^ / cos(k9 m ■ (x c - x s ))J (k\x c - x s \)dk 



1 



iV ^ I , -> \ 2 

m=1 I- -1/1 / Xc ~~ Xs 

X f X q 



(9) 

Theorem 3.3 (Case 2: Non-symmetric incident directions). Under the configuration 
of a small number of non-symmetric incident directions, the structure of fifjj) becomes 



i • —X ( n Xc Xs 

M 1 sm I m ■ — — 

y X c X s 



m=1 l" I la Xc ~ Xs 



\X C X s 



for a sufficiently large F and kp — > oo. 



Proof. In constrast with Theorem I3.2[ we assume that the incident directions are non- 
symmetric. Then, (J7J) can be written as 



^td(^ s ; k f ) 

=3? 



. M 



— V" exp(ik f 9 m ■ (x c - x s ))$(x c , x s ; kf) 

J m=l 
' 1 M 

— ^2 exp(ik f e m ■ (x c - x s )){J (k f \x c - x s \) + iY (k f \x c - x s \)} 



f 

J m=l 
M 

cos(k f 6 m ■ (x c - x s ))J (k f \x c - x s \) 



m=l 



- sxxx(kfO m ■ (x c - x s ))Y (k f \x c - x s \) 



and therefore 



^td{x s ] kf) 



M 



xnzx\R TD (x s -k f )\ 2M^ 



cos(k f 9 m ■ (x c - x s ))J (kf\x c - x s \) 
- sxn(k f 9 m ■ (x c - x s ))Y (k f \x c - x s \) 
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Let us recall an integral formula in 0, formula 11.4.40 (p. 487)], 

2i 



exp(iat)Y (bt)dt = < 



2i 



In 



|2 /~r2 



sin 



-i 



Va 2 - b 2 



K Tv\Za 2 — b 2 
Then, if kp — > oo, we can evaluate 

Y] sin(A; f 6 m ■ (x c - x s ))Y (k f \x c - x s \) « / 
/=i ^ 



1 



Va 2 - 6 2 



if < a < b 



if < b < a. 



sm(k9 m ■ (x c - x s ))Y (k\x c - x s \)dk 



- sin 1 \9 m - —; 



Q / exp(k6 m ■ (x c - x s ))Y (k\x c - x s \)dk 



71 



0C c 'X a 



0C c 0C q 



(10) 



\x c - x s \\\l- ( 9, n 



x c x s 



Hence, by combining dH), ©, an d (|10|) . we can obtain 



MVITDl^s 



{x s ; k f ) 



1 M 

F ^ max |R TD (f s ; fc ; ) | ~ 2M ^ 



1 - - sin" 1 I ft 



7T 



3j C Xi C 



Xn X m q \ \ / 1 ( 9m, ' 



\x c x s \ 



(11) 

Theorem 3.4 (Case 3: Sufficiently large number of incident directions). Under the 
configuration of a sufficiently large number of incident directions, the structure of |Qj be- 
comes 

A(\x c — x s \; kp) — A(|5? c — x s \; ki) if kp < +oo, 
(5(x c , x s ) 

for an adequately large F. Here, A(x) is defined as 

A(x; k) := k( J (kx) 2 + J^kxf 



^MTD 



%; F) 



if kp — > oo, 



(12) 

and 5(x c ,x s ) denotes the Dirac delta function. 

Proof. Suppose that M is large enough, i.e., M — > oo, then by [lil Lemma 4.1], ((?]) can 
be approximated as 



°td(£ s ; k) ~ 3? 



3? 



7T 



^ exp(zA;(9 m • (x c - x s ))®(x c , x s \ k) 



m=l 



2k 



— yJ (k\x c - x s \) + iY (k\x c - x s \) ) I exp(ik6 ■ (x c - x s ))d0 
J (k\x c - x s \) 2 . 



Assume that F is large enough and that Uf < +00. Then applying an integral formula 



of the Bessel function [36|, page 35] 

J {t) 2 dt = t[ J {t) 2 + J^tf ) + / J x {tfdt 



and a change of variable t = k\x c — x s \ yields 

F -y i>k F j rk F \x c -x s \ 

Vl TD (^; k f ) ~ t ^ / Jo(kWc ~ x s \fdk = — — / J (t) 2 dt 

, -1 ^F — J^lJk, \ x c ~ x s I Jk-I \Sr.SJ 



\2. 



/=1 

=A(|x c - x s |; fci?) - A(|x c - x a |;fci) + / Ji(£;|f c - £ s |) 2 d/c, 



(13) 



where A(t; fc) is given in ffl2]) . Note that the term 



F 



6(|x c — k) := / Ji(fc|x c — x a \) dk 



can be negligible since A(\x s — x c \, k F ) = 0(k F ) and Q(\x c — x s \, k) <C 0(kp), refer to [20 



Finally, by applying kp — > 00 in f TT3"j) . we can immediately obtain the following result: 
f> TD (f s ; A;,) ~ J™ J (k\x c - x s \fdk = I + 00 |[ % = J (14) 

Remark 3.5. Based on structures fifty, / T77]) . and / T73j) . we can summarize some properties 
of (0|) as follows: 

(PI). Formulas fi9\) and 01]) reaca iaezr maximum value at x s G Z> iaat satisfies 

I *^ C *^ s I — * 

Therefore, at x s = x c G C, tae map 0/ (0|) yieid 1 iaai tae location of C can be 
accurately determined. However, at point x s G T> such that 

\x c x s 1 

iae map of fifty will show a large magnitude such that unexpected replicas C will appear 
(see Figure^). 
(P2). If x s satisfies 



~1, 



@m ' \ x c x sJ 

then both the numerator and denominator of fiTI\) become 0. Hence, at this point, fiTI\) 
will yield a large value at not only x c but also other locations in T>. Therefore, the 
configuration of symmetric incident directions will yield a better imaging result (see 
Figure^). 

8 



(P3) 



If the number of incident directions is increased to as many as possible, i.e. 



M 



oo 



it is clear that a more accurate location of x c e C can be obtained. In this case, an odd 
number of M will not affect the imaging performance significantly. 
(P4). Application of multiple frequencies enhances the imaging performance. 



Figure 1: Illustration of Remark 13.51 At the location x s such that x s = x c (black colored vector) and 
9 m ■ i| c ~| s i = ±1 (green colored vector), true shape of C and ghost replicas C will appear, respectively. 



4. Numerical experiments 

In this section, we exhibit some numerical examples to demonstrate the effectiveness of 
the proposed algorithm. Throughout this section, the homogeneous domain D is chosen as 
the interior region of the two-dimensional unit disk centered at the origin, i.e., 



The adopted wavenumber has the form kf — ^ for / = 1, 2, • • • , F; here, A/ is the given 
wavelength that is equi-distributed between Ai = 0.7 and A^ = 0.4. Based on Theorem 13.21 
the incident direction 6 m is selected as 




O 



B = {(cost, sin t) : t G [0,2tt]}. 





we choose two single cracks C\ and C 2 as 



C 2 



c 



1 — 




<*±^): t6 [-l,l]} 



t e [-1 
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and inspired by [34|, we choose multiple cracks Cm as 
C M = {(t - 0.2, -0.5t 2 + 0.5) : t E [-0.5, 0.5]} U {(t + 0.2, t 3 + t 2 - 0.5) : s E [-0.5, 0.5]}. 

Note that in order to show the robustness of the proposed algorithm, a white Gaussian 
noise with 15 dB signal-to-noise ratio (SNR) is added to the unperturbed data. 

First, let us consider the imaging result of C\. The left-hand side of Figure [2] shows the 
map of IRmtd (5^ ; 4) for M — 16 different directions. In this result, because of its small 
magnitude at (0.46,0.42), C\ appears divided. Moreover, owing to the large magnitude 
at certain points (for example, (0.23,0.18)), we cannot recognize its shape at this stage. 
The right-hand side of Figure [2] is the map of Rmtd(^; 16). The figure shows that there is 
undivided or whole C\\ this means that we obtained a more accurate result than the previous 
one. However, because certain points have a large magnitude, we cannot determine the true 
shape of C\. 
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Figure 2: Map of R M tdO? s ; F) with F = 16, M = 4 (left) and F = 16, M = 16 (right) for C x . 



Figure [3] shows the map of Mmtd(^s! F) for crack C<i- In contrast with the result of Ci, 
it is difficult to identify the shape of the crack because of the great number of points with 
large magnitudes. Fortunately, similar to the results in Figure EJ map of Mmtd(^ s ; 16) yields 
a better image than that obtained by the map of Mmtd^s; 4). 

The application of multiple cracks Cm is shown in Figure HI This figure shows the maps 
of RtdO? s ; fij) and Mmtd(^ s ; F). Throughout these results, we can observe that the map 
of Mtd(^ ? s; k) yields a poor result even if the number of incident directions M is large; 
conversely, the multi-frequency imaging function Mmtd^; F) yields an accurate result even 
if M is small. 

In FigureEl maps of Mmtd(^s! 16) with M = 5 number of incident directions are displayed 

for cracks C\ and Cm in order to verify the (P2) of Remark 13.51 By comparing Figures |5] and 

HI it is difficult to discern the shape of true crack(s) because of the appearance of unforseen 

replicas with large magnitudes. Note that if M is an odd number but is sufficiently large 

enough, the map of 1mtd(^ s ; F) yield a good result, refer to Figure [6j 
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At this moment, we increase the number of applied frequencies F and incident directions 
M. Figures [7J [BJ and |9] depict the map of Kmtd(^ s ; F) when F = 64 and M = 16 or 



M = 64. As we expected in (P3) of Remark 1331 the results yielded by Kmtd(i s ; F) improve 



with increasing number of F and/or M. 

Remark 4.1 (Limitation of the multi-frequency topological derivative). Although multi- 
frequency topological derivative imaging yields very good imaging results for not only single 
but also multiple cracks, it still has the following limitations. 

1. Let us consider the imaging of a crack of a large curvature near one of the crack tips. 
To illustrate this, we choose a crack from \23t] : 

C 3 = {(t + t 2 ,0.5t 4 + t) : -1 < t < 3} . 

From the results in FigureXUX the imaging results yielded by M MTD (x s ; F) are satisfac- 
tory. However, at the point of the large curvature, the map o/1mtd(^s! F) reconstructs 
a slightly different shape of the crack. 

2. We apply the proposed algorithm for imaging an oscillation crack represented as fol- 
lows: 

C 4 = { (t, 0.5t 2 + 0.1 sin(37r(t + 0.7))) : -0.6 < t < 0.6} . 

FigureUl\ shows the map o/Rmtd(^s! F) . The obtained result by the proposed algorithm 
does not improve even when the number of F and M are increased. The proposed 
algorithm requires further optiization. 



5. Concluding remarks 

We improved the traditional topological derivative concept to determine location and 

reconstruct the shape of arbitrary-shaped curve-like perfectly conducting cracks. For this 

11 



x-axis 



x-axis 




purpose, we designed a reconstruction algorithm based on the multi-frequency topological 
derivative, analyzed its structure, and investigated its properties. We confirmed that the 
multi-frequency topological derivative contains singularity on the cracks, and therefore it 
reaches its maximum value at the location of the cracks. Moreover, we demonstrated that 
why ghost replicas appeared in the neighborhood of cracks and why the configuration of 
symmetric incident directions yields better results than the non-symmetric configuration. 
Numerical simulations under various situations depicted both the benefits and limitations of 
the proposed algorithm. Because the algorithm requires only a one-step iteration procedure, 
it is very fast in terms of shape identification but does not yield the complete shape of 
cracks. However, by adopting the obtained result as a starting point of a Newton-type 



based algorithm [22J, it is expected that complete shape reconstruction can be successfully 
accomplished. 

This paper deals with the reconstruction of perfectly conducting cracks with Dirichlet 
boundary condition. Accordingly, the extension of this research to the crack with Neumann 
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Figure 5: Map of Rmtd^s;^) with F 
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Figure 6: Map of R M td(^; F) with F 




-1 -0.5 0.5 1 



x-axis 

16, M = 5 for d (left) and C M (right). 
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16, M = 15 for Ci (left) and C M (right). 



boundary condition will be an interesting work. We believe that applying higher-order terms 



in the asymptotic expansion formula [i| give a higher order topological derivative [13| . The 
calculation and analysis of higher-order topological derivative will be a valuable research 
topic. According to the Statistical Hypothesis Testing jf|, multi- frequencies are expected to 
enhance the imaging performance. The corresponding work will include a careful stability 
and resolution analysis of the multi-frequency topological derivative. 
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Figure 10: Map oi R M tb(x s ; F) for C 3 with F = 16, M = 16 (top, left), F = 16, M = 64 (top, right), 
F = 64, M = 64 (bottom, left) and true shape (bottom, right). 
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Figure 11: Same as Figure HU1 except the crack is C4. 



[35] W.-K. Park and D. Lesselier, Reconstruction of thin electromagnetic inclusions by a level set method, 

Inverse ProbL, 25 (2009), 085010. 
[36] W. Rosenheinrich, Tables of Some Indefinite Integrals of Bessel Functions, 2011, available at 

http : //www. fh-jena.de/~rsh/Forschung/Stoer/besint .pdf 
[37] B. Scholz, Towards virtual electrical breast biopsy: space frequency MUSIC for trans-admittance data, 

IEEE. Trans. Med. Imaging, 21 (2002), 588-595. 
[38] J. Sokolowski and A. Zochowski, On the topological derivative in shape optimization, SIAM J. Control 

Optim., 37 (1999), 1251-1272. 



17 



